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We demonstrate the tip induced control of the spin state of copper phthalocyanine (CuPc) on an 
insulator coated substrate. Accounting for electronic correlations, we find that, under the condition 
of energetic proximity of neutral excited states to the anionic groundstate, the system can undergo a 
population inversion towards these excited states. The resulting state of the system is accompanied 
by a change in the total spin quantum number. Experimental signatures of the crossover are 
the appearance of additional nodal planes in the topographical STM images as well as a strong 
suppression of the current near the center of the molecule. The robustness of the effect against 
moderate charge conserving relaxation processes has also been tested. 

PACS numbers: 85.65.+h,68.37.Ef,73.63.-b,75.30.Wx 


Introduction - Research on single molecule junction- 
shas witnessed in recent years a broadening interdisci¬ 
plinary interest pQ. For example, spin dependent trans¬ 
port mm or nuclear spin resonance [4] have been inves¬ 
tigated. In this emergent field of molecular spintronics, 
spin-crossover metalorganic compounds (SCOs) play a 
prominent role nm. These molecules undergo a transi¬ 
tion between metastable spin states under the influence 
of external stimuli m ■ The many-body exchange inter¬ 
action of the ^-electrons on the metal center, in combi¬ 
nation with the crystal field generated by the surround¬ 
ing ligand, determines their spin state. In three-terminal 
devices, the change of charge state tuned by the gate 
electrode has been shown to govern the associated spin 
state Recently, SCOs have been in the focus of 

STM experiments mmm- More generally, the role of 
many-body effects in STM single molecule junctions is 
receiving increasing attention, both theoretically [IMS] 
and experimentally MEEIEE7]. 

In this Letter we demonstrate the appearance of a non¬ 
equilibrium high-spin state in CuPc on an insulating sub¬ 
strate caused by population inversion, and show experi¬ 
mentally observable fingerprints of this effect. We illus¬ 
trate that, for a given substrate work function, it is pos¬ 
sible to control the effective ground state of the molecule 
by varying the tip position or the bias voltage across the 
junction. The only requirements for this genuine many- 
body effect are an asymmetry between tip and substrate 
tunneling rates, which is naturally inherent to STM se¬ 
tups, and an energetic proximity of an excited neutral 
state of the molecule to its anionic ground state. As dis¬ 
cussed below, the experimental set-up is similar to that of 
Ref. [17], but with a slightly larger workfunction for the 
substrate. Control over the workfunction can be achieved 
by choosing different materials or crystallographic orien¬ 
tation for the substrate, with effects analogous to a dis¬ 
crete gating of the molecule. Several approaches to gate 
an STM junction have been also very recently investi¬ 
gated [18U2Q1 . 

Many-body Hamiltonian and spectrum of CuPc - To 
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FIG. 1. (Color online) (a) Frontier orbitals used for the many- 
body calculation, in their complex representation. The color 
code shows the phase of the wavefunctions. (b), (c) Full and 
low-energy cutout, respectively, of the many-body spectrum 
of CuPc at chemical potential p — —4.65 eV. (d) Scheme of 
the lowest-lying many body states. 


properly describe the many-body electronic structure of 
CuPc is by itself a nontrivial task, since the relatively 
large size of the molecule makes it impossible to diagonal¬ 
ize exactly a many-body Hamiltonian written in a local, 
atomic basis as done for smaller molecules [2TH23] . STM 
transport experiments on single molecules, however, are 
restricted to an energy window involving only the low- 
lying states of the molecule in its neutral, cationic and an- 
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ionic configuration, with the equilibrium configuration at 
zero bias set by the workfunction 0 O of the substrate El- 
This allows one to use a restricted basis of frontier or¬ 
bitals to construct the many-body Hamiltonian [24] . For 
example, for a copper substrate as in [17] is 0o = 4.65 eV, 
and CuPc in equilibrium is in its neutral ground state. 
Thus, in the following we only retain four frontier orbitals 
of CuPc, the SOMO (5), the HOMO (H) and the two 
degenerate LUMO (L^) orbitals, see Fig. 1(a). In equi¬ 
librium, the molecule contains Nq = 3 frontier electrons. 
In this basis, all matrix elements of the Coulomb interac¬ 
tion are retained. Hence, besides Hubbard-like density- 
density interaction terms, our model also includes ex¬ 
change and pair hopping terms, which ultimately are im¬ 
portant for the structure and spin configuration of the 
molecular excited states. 

The Hamiltonian of CuPc in the basis of the four single 
particle frontier orbitals reads 

Hmol = £i hi + ^ Vijkl ^ia^ka'^la'^-ja ^ (1) 

i ijkl a a' 

where i = 5, iJ, L± and a is the spin degree of free¬ 
dom. The energies e* = e* + contain the single parti¬ 
cle molecular energies e* obtained from diagonalizing the 
single particle Hamiltonian Ho of CuPc, es = —12.0 eV, 
€h = —11.7 eV and €l± = —10.7 eV. The parameters A i 
account for crystal field corrections and the ionic back¬ 
ground of the molecule, since the atomic onsite energies 
in H 0 come from Hartree-Fock calculations for isolated 
atoms [25] . The A i are free parameters of the theory. 
Isolated CuPc has D^ symmetry; the four molecular 
orbitals | ia) that make up the basis of Eq. 0 trans¬ 
form like its big (5), a\ u ( H ) and e u (T ± ) representa¬ 
tions. As a consequence, they acquire distinct phases 
0i when rotated by 90 degrees around the main sym¬ 
metry axis of the molecule, as illustrated in Fig. [lja). 
This yields an easy rule to determine the nonvanishing 
Coulomb matrix elements Vijki in Eq. 0 : V ijkl ± 0 if 
0i — (j>j + (\>k ~ <\>i — 0 mod 27r, i.e., nonvanishing contri¬ 
butions are only possible if the phases of the correspond¬ 
ing molecular orbitals add up to multiples of 27r. These 
considerations remain true in the presence of a homoge¬ 
nous substrate which reduces the symmetry to C± v . For 
a detailed discussion concerning the parametrization of 
Eq. 0 we refer to the supplemental material [26] . Ex¬ 
act numerical diagonalization of H mo i finally yields the 
many body eigenenergies E^m and eigenstates | Nm) of 
the molecule, labelled after particle number N and state 
index m. 

Since the molecule is in contact with the substrate and 
is able to exchange electrons, it is necessary to consider 
a grandcanonical ensemble H mo i — /UV, where (i is the 
chemical potential of the substrate which is given by its 
negative workfunction, p = —0o- Moreover, the pres¬ 
ence of the leads renormalizes the Hamiltonian Hq due 


to image charges effects nanz]. We model these effects 
with an effective Hamiltonian H mo i-env = — 5{ C (N — Aq) 2 , 
with N the particle number operator on the system and 
S[ c obtained from electrostatic considerations, see supple¬ 
mental material. To fit our spectrum to the experiment 
of Swart et al. m , which was taken on a copper substrate 
Cu[100] (0o = 4.65 eV) on a trilayer of NaCl, we used a 
constant shift A^ = A = 1.83 eV, a dielectric constant 
e mo i = 2.2 in the evaluation of the matrix elements 
and an image-charge renormalization S- 1C = 0.32 eV. 

Figures[ljb), (c) show the cationic, neutral and anionic 
subblocks of the many particle spectrum and their degen¬ 
eracies. The neutral groundstate has a doublet structure 
(with total spin S = |) coming from the doubly filled 
HOMO and the unpaired spin in the SOMO, whereas 
the cationic and anionic groundstates have triplet struc¬ 
tures (5 = 1). The former has a singly filled HOMO, 
the latter a singly filled LUMO orbital which form spin 
triplets (and singlets, 5 = 0, for the first excited states) 
with the singly filled SOMO. Finally, the orbital degen¬ 
eracy of the LUMO makes up for an additional twofold 
multiplicity of the anionic ground and first excited states. 
The first excited state of the neutral molecule is found to 
be also a doublet (5 = with additional twofold orbital 
degeneracy. Finally, the second excited state shows a 
spin quadruplet structure (5 = §) together with twofold 
orbital degeneracy. A schematic depiction of these states 
is shown in Fig.Qd). As the actual states are linear com¬ 
binations of several Slater determinants, only dominant 
contributions are shown. 

Transport dynamics and spin-crossover - The full sys¬ 
tem is characterized by the Hamiltonian H = H mo i + 
Hmoi-env + H s + H T + H tun , where H s and H T are 
describing noninteracting electronic reservoirs for sub¬ 
strate (S) and tip (T). The tunneling Hamiltonian is 

Htun = *ki 4kAa + hx -’ where creates an 

electron in lead p with spin a and momentum k. The 
tunneling matrix elements tjL are obtained analogously 
to Ref. [12]. The dynamics is calculated via the Gener¬ 
alized Master Equation for the reduced density operator 
p re d = Trs,T (p), see Refs. [12, 22]. In particular, we are 
interested in the state p^ d solving the stationary equa¬ 
tion jC [p r ed] = 0, where C is the Liouvillian superopera¬ 
tor. 

In analogy to Ref. [28] , we included a phenomenological 
relaxation term C ve \ in the Liouvillian [29] : 

£rel [p\ = “I ip -^Pmm \ Nm ) ( Nm \ E Pnn ) • ( 2 ) 

\ Nrri n / 

It is proportional to the deviation of the reduced den¬ 
sity matrix from the thermal one p th , which is given by 

the Boltzmann distribution p^m ~ ex P "ffr ) w ^h 

Pmm = 1- Since C Te \ describes relaxation processes 
which conserve the particle number on the molecule, it 
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cationic resonance: 0o = 4.65 eV 
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FIG. 2. (Color online) Constant height current maps (a,d,g), 
constant current maps (b,e,h) and maps of the system’s to¬ 
tal spin S (c,f,i). Constant height and spin maps are taken 
at a tip-molecule distance of 5 A, constant current maps at 
currents I — 0.5, 0.75, 1.0 pA for panels (c), (f), and (i), re¬ 
spectively. 



does not contribute directly to the current. The relax¬ 
ation factor 4 is taken of the same order of magnitude of 
the tip tunneling rate. The stationary current through 
the system is evaluated as 

(Is + It) = ^ (N) = Ttmoi (VC^]) = 0, (3) 

The Liouvillian C = £ re i + decomposes into the 

relaxation term and sub-Liouvillians for each lead. Sort¬ 
ing of the occuring terms in Eq. © after substrate and 
tip contributions yields the current operator of the re¬ 
spective lead r] as I v = NC ?7 . 

Results of our transport calculations are presented in 
Fig. [2j In panels (a,d,g) we show constant height cur¬ 
rent maps, constant current STM images in (b,e,h) and 
in (c,f,i) maps of the expectation value of the total spin 
of the molecule depending on the tip position, S rT = 

\J (£ 2 )r T + i - | where ( S 2 ) rT = Tr mol (5 2 ^(r T )) • 
The constant height and spin maps are each taken at 
a tip-molecule distance of 5 A. The upper three pan¬ 
els (a,b,c) are for a workfunction of 0o = 4.65 eV and 
a bias voltage of V& = —2.72 V. At this position the 
cationic resonance is occuring. Since the difference be¬ 
tween neutral and cationic groundstate is the occupation 
of the HOMO (see Fig. [j|d)), tunneling occurs via this 
orbital and the current maps (a,b) resemble its structure. 
With the same work function 0o = 4.65 eV, the anionic 
resonance is taking place at positive bias Vb = 0.81 V, 
see Fig. ©d, e). For equivalent reasons as in the former 
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FIG. 3. (Color online) (a) Differential conductance and (b) 
total spin curves taken at different tip positions and work- 
functions around the bias V r es(0o) of the anionic resonance. 
The inset in (b) shows the change of the spin for the standard 
case in magnification, (c) Populations of the density matrix 
around V r es(0o)- Left panel: standard case, 0o = 4.65 eV. 
Middle (right) panel: anomalous case, 0o = 5.2 eV, with tip 
near the center (outer on the ligand). 


case, tunneling is happening via the LUMO and the spa¬ 
tial dependence of the current resembles the topography 
of this orbital. Panels (g,h,i) are recorded, instead, at 
0o = 5.2 eV, again at the anionic resonance which is 
now shifted to Vb = 1.74 V due to the larger workfunc¬ 
tion. Panel (g) is puzzling. Despite being an anionic 
resonance, it closely resembles the HOMO, cf. as pan¬ 
els (a)-(b). A closer inspection, though, reveals also an 
alikeness with the LUMO (see panel (d)) but with addi¬ 
tional diagonal nodal planes, matching the nodal plane 
structure of the HOMO. When observing in panel (h) 
the constant current map, and comparing it with panels 

(b) and (e), this statement becomes more evident. This 
anomalous topography can not be explained by single 
orbital tunneling. 

Panels (c), (f) and (i) reveal the tip-position depen¬ 
dent expectation value of the total spin. At the stan¬ 
dard anionic transition, panel (f), the spin remains essen¬ 
tially constant. At the standard cationic transition, panel 

(c) , a rather homogeneous enhancement of the molecular 
spin is due to small populations of a large number of ex¬ 
cited states, made accessible by the large resonance bias 
(V res = —2.7 V). The anomalous anionic transition, panel 
(i), shows, however, the largest variation of the molecu¬ 
lar spin, concentrated at the positions of the anomalous 
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FIG. 4. (Color online) Simplified sketch of the tunnelling 
processes at the anionic resonance for the standard (0o = 
4.65 eV) and the anomalous (0o = 5.2 eV) case. In the latter 
population inversion takes place. The colors of the arrows 
denote tip positions where the corresponding transition acts 
as a bottleneck: Orange (blue) stands for the center (the outer 
ligand) of CuPc. 


current suppression, compare panels (g) and (d). To ex¬ 
plain the unconventional properties shown in Fig. [2j we 
examine bias traces taken at different tip positions and 
values of the workfunction. Figure J3fa) shows a shift of 
the anionic resonant peak in the for the anomalous 
case. The value V res at which the peak is expected is 
given by 

Kes(0o) = -j—f (^C/Vo + 1,0 — En o ,0 ~ ^ic + 00 ) 5 (4) 

a T \e\ 

where is the fraction of bias drop between tip and 
molecule, and E ^,o is the energy of the 7V-particle ground 
state. The shift of the resonance to lower biases seen in 
Fig-i a) suggests the appearance of a population inver¬ 
sion from the neutral ground state to an excited state. 
Transitions from the latter to the anionic ground state 
open in fact at much lower biases. Also the evolution 
of the spin of the molecule shown in Fig. [3jb) reinforces 
this proposition. In the anomalous case, the change of the 
system from a low to a high spin state, as well as the satu¬ 
ration of the spin, can be clearly seen. This contrasts the 
normal anionic transition, where only a marginal change 
is observable. In Fig. [3jc) we show the evolution of the 
eigenvalues of the stationary density matrix p^ d , i.e. the 
populations of the physical basis |22], around the anionic 
resonance V^ es (0o)? depending on workfunction and tip 
position. In the standard case (left panel of Fig. He)), the 
ground state of the system is always the neutral ground 
state. For the anomalous case (middle and right panels 
of Fig. [3| however, the picture changes dramatically, as 
there is a remarkable depopulation of the neutral ground 
state in favor of different excited states, depending on 
the position of the tip. 

We focus now on the mechanism explaining population 
inversion and associated spin-crossover. In the standard 
case, at sufficiently high bias, the transition from the neu¬ 
tral to the anionic groundstate is opening, and tunneling 
of an electron into the LUMO brings the molecule into 


the anionic ground state. By consecutive tunneling to the 
substrate, the system goes back into its neutral ground 
state, see Fig. [4] for a simple sketch. Since the tunnel¬ 
ing rates to the substrate are much larger than their tip 
counterparts, the system stays essentially in the neutral 
ground state with spin S «= \. Also in the anomalous case 
an initial tunneling event brings the molecule into the 
anionic ground state. However, from there, due to finite 
temperature and proximity of the many-body eigenen- 
ergies, the system has a finite probability to go into a 
neutral excited state by releasing an electron to the sub¬ 
strate. The position of the tip and the structure of these 
excited states themselves then determine the stationary 
state: The molecule can only return to its neutral ground 
state by successive transitions to the anionic ground state 
via the tip, and from there to the neutral ground state 
via the substrate. However, the former process acts as 
a bottleneck and depends on the tip position. Leaving 
the first excited state (S = ^) requires tunneling into the 
SOMO, while leaving the second excited state (S = |) 
would require tunneling into the HOMO. Additionally, 
near the center of the molecule the HOMO is vanishing, 
whereas on the outer ligand part the SOMO has little to 
no amplitude. Therefore, tunneling into these orbitals at 
the respective positions is strongly suppressed and the 
system ultimately ends up in the corresponding neutral 
excited states. 

Conclusions - For an experimentally accessible sub¬ 
strate workfunction of 0o = 5.2 eV, we predict the ap¬ 
pearance, in proximity to the anionic resonance, of pop¬ 
ulation inversion between the neutral ground and ex¬ 
cited states of CuPc. Depending on the tip position, the 
molecule is triggered into a low-spin (S=l/2) to high-spin 
(S=3/2) transition which is mediated by this population 
inversion. This inversion is experimentally observable 
via dramatic changes in the topographical properties of 
constant height and constant current STM images, com¬ 
pared to a standard LUMO-mediated anionic transition. 
Direct observation of the spin-crossover might be accessi¬ 
ble using spin-polarized scanning probe microscopy tech¬ 
niques. m The effect is also robust against moderate 
charge conserving relaxation processes. The quantita¬ 
tive accuracy of the spectroscopic and topographical re¬ 
sults presented in this Letter is limited by the adopted 
semiempirical model. The spin-crossover with the asso¬ 
ciated anomalous topography of the anionic resonance 
depends, however, on qualitative properties of the many- 
body spectrum and of the molecular orbitals. Thus, de¬ 
spite our focus on CuPc, they should be observable also 
in other molecules with comparable frontier orbital struc¬ 
ture. 

The authors thank Thomas Niehaus, Jascha Repp and 
Dmitry Ryndyk for fruitful discussions. Financial sup¬ 
port by the Deutsche Forschungsgemeinschaft within the 
research program SFB 689 is acknowledged. 
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THE SINGLE PARTICLE HAMILTONIAN OF CUPC 


The single particle Hamiltonian of CuPc without many-body interaction can be written in second quantization as 


Ho ^ ^ ^ara,/3n d arn(7 dp n(J 

a/3cr 

rnn 


(i) 


where d arna creates an electron in the atomic orbital \ama) with quantum numbers m, a centered at atom a: 

(r| ama) = ( 2 ) 

The matrix elements are defined as 

ham,(3n • = ^am^a/3^mn T b(xm,f3m (3) 

where e arn is the energy of orbital m on atom a , and b arn ^ n is the hopping integral between orbital m on atom a and 
orbital n on atom (3 . The atomic onsite energies are taken from Ref. [1] and the hopping integrals are calculated in 
nearest-neighbour approximation following the LCAO schemes proposed by Refs. [2, 3], analogously [4] to our earlier 
work on H 2 PC [5]. Geometrical parameters needed for the LCAO calculation like interatomic distances and angles 
are taken from Ref. [6]. For the ligand we consider the set of all 2s (Is for hydrogen), 2pa, and 2p y orbitals as the 
cr-system, and consequently the set of 2p^ orbitals as the 7r-system. On the metal the 3d xy , 3 d x 2 _ y 2 , 3 d z 2 and 4s 
orbitals contribute to the cr-system, while the ?>d zx and 3d yz belong to the 7r-system. The single particle energies e* 
and molecular orbitals \ia) of the CuPc molecule are obtained by numerical diagonalization of Hq- 


SETTING UP THE MANY PARTICLE HAMILTONIAN 


The full many-body Hamiltonian of a molecule in Born-Oppenheimer approximation in second quantization reads 
H-mol Hq T Vee ^ ^ ^ara,/3n d~ ~ ^ ^ ^ ^aB^yS ^ama^'yp<j'^-5q<j'^-/3n<j , > (^) 


a/3a 

rnn 


cx.fB ^/5 crcr' 
mnop 


where V™p™° s p is the matrix element of the Coulomb interaction, 

V afh°8 P = J d 3r l J d 3 f 2 C( r l) 4>Pmr(ri) ^ | <^oa' ( r 2 ) 4>Sp*> (r 2 ). (5) 

After performing a transformation to the molecular orbital basis, in which Hq is diagonal, d| a = J2am Cam dl mcr , and 
using the approximation that the basis |<amcr) is orthogonal, i.e. ( amcr\f3ncr') = ^a/^mn^crcr 7 ? Eq. (4) reads: 


H mo l — ^ ^ ^icr d~ EE Vijki dj, 


dL'd,„,d 


icr U /ccr / U ^cr / U jcr 5 


ijkl crcr' 


(6) 


where the matrix element of the Coulomb interaction now reads 

Vijki = j d3 G J d 3 r 2 ip*(ri)tpj(ri) ^ ^(r 2 )^i(r 2 ) ■ 


(7) 


To cut down the size of our Hilbert space we split the basis into frozen and dynamic states, where Nf frozen states 
are assumed to be always fully occupied and N e set to be always empty, whereas we do not make any assumption 
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about the occupation of the Nj dynamic states. In occupation number representation a general state of the Fock 
space then looks like 


\W) « 1 11. . . 11 n k \n ki ... n; T n 4 pO. . . 00 ). 

2N f " 2 N, / 2 N e 


(8) 


The different terms of the Coulomb interaction 


The Coulomb interaction 

( 9 ) 

ijkl < 7 <j' 


can be split into different terms corresponding to the number of different indices in the sum. The first term, where 
all four indices are the same is: 


^ ^ Van 


( 10 ) 


Terms with two different indices are: 


2 Viijj fliflj + (Viiij flier di^djp + h.C.^ 2 Vijji 

[ij] [ij]a [ij]cr 


flier fljcr ^ia^ja^ia^-ja 


) + )£ 


Vijij ^icr^ia^ja^jcn 


(ii) 


where we used that Vijki = VkUj — Vj*uk = ^ikji • The symbol [...] under the sum means for example for [ijk\: 
i 7^ j,i 7^ k,j ^ k. Terms with three different indices are: 


[tjk] 


ijik ^ia^ia^ka^j a + h.c.^ + E^ E^ Vijki ( 

[ijk] cr 


^ia^ka^ia^ja ^lia ^ka^ja 


) + EE 

[ijk] cr 


Vujkhi d] a d k(7 . (12) 


Consequently the last term with four different indices reads 


2 ^ ^ > Vijki ^jg^kcr'^la'^ja’ 

[ijkl] era' 


(13) 


Symmetry considerations 


The molecular orbitals | i) are transforming under symmetry operations of the group of the molecule according to 
its irreducible representations and are thus acquiring distinct phase factors. Let R be a rotation of 90 degrees around 
the fourfold molecular axis with R~ 1 r =: r: 


R Ij) = e*' | j ), 

V’j(^ _1 r) = tpj( f) = (r|j) = (r|i?|j) = (r| j) = ipj( r)e*^. 


(14) 


This yields for Vijki after a coordinate transformation r i under the integral: 


Vijki = J d3ri J d3r 2^*(ri)^j( r i) | ri Z ~^ k( r 2)M r 2) 

= inei J ^ / d3f2 ^( ?1 )^( ?1 ) | fl Z r 2 | 

= J d3ri / d3r 2 'i’i {R~ lr i)^Pj(R~ lr i) | r ^ ^ r2 | ^fc(r 2 )^(r 2 ) 

= ^~ f d 3 n f d 3 r 2 V’*(ri)e-^^ i (r 1 ) e ^ 1 --- 

47reo J J |ri-r 2 | 

= e -i(<Pi-<Pj+<Pk-<Pl) y. jkh 
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Us 

11.352 eV = 

— 7 P 

J H-\ — 

548 meV 

U H 

1.752 eV 4 X _ 


258 meV 

1 

II 

1.808 eV J p _ 


168 meV 

USH 

1.777 eV J e s l = 

1 

+ 

1 

9 meV 

U SL 

1.993 eV J e s *H = 

J P SH 

2 meV 

Uhl 

1.758 eV 




TABLE I. Major nonvanishing Coulomb integrals between the SOMO(S'), the HOMO (if), the LUMO + (+) and the LUMO - 
(-). All values are calculated numerically using Monte Carlo integration [7] of the real space orbitals depicted in Fig. 1(a) of 
the main publication and renormalized by a constant e mo i = 2.2. 


In order to obtain a true statement the phases have to meet the constraint 

fa - fa + (j) k - fa = 0 mod 27r. (15) 

This condition drastically reduces the number of possible matrix elements of the Coulomb interaction, since the SOMO 
(S), the HOMO (H) and the two degenerate LUMO (=b) orbitals gather different phases under rotations of 90 degrees 
around the main symmetry axis of the CuPc molecule: 


<t>S = 7T, 

4>h = o, 

^=4 

This we will exploit to further simplify Eq. ((6)). Thus, for Vijki it holds that: 

V ijk i = v ijkl 

=> Vijki 4 o if fa - 4>j + (j)k - 0 mod 27r, 

since all phases fa are different; fa 4 4>j f° r ^ 4 T Putting all together, Eq. ((6)) can be recast in the form 


Hmol — ^ ^ (C ^i) fli T ^ ^ Ui Tli^Tli^ T ^ ^ ^ Vij fliflj ^ EE 

i i [ij] [ij] a 


_ jt jt j ^ 


+ 2 E E J ij btdLdjAv + 2 E E ( 

[ij] a [ijk] a 


[ij] a 

^ijk ^ia^ia^ka^j a "f~ ^ 


) 


-) + )£E Vijki 


ia^-ka' ^ la' ^jcr i 


[ijkl] crcr' 


(16) 

(17) 

(18) 


(19) 

( 20 ) 


( 21 ) 


where indices are now running only over orbitals from the dynamic set (S', iL, L + , L~). The abbreviations we 
introduced in Eq. ((21)) are the orbital Coulomb interaction Ui = Van, the inter-orbital Coulomb interaction Uij = 
Vujj , the exchange integral Jf* = Vijji , the ordinary pair hopping term JP. = Vi jij and the split pair hopping 
term j?. fe = V^/c- Matrix elements 4j/cZ of the Coulomb interaction are calculated numerically by Monte Carlo 

integrating [7] the respective single particle orbitals with the Coulomb interaction V^iT,^) = (47 T£o |it — r* 2 1) _ 1 and 
renormalizing the bare values by a constant screening factor e mo \ = 2.2 to account for the screening given by the 
electrons in the frozen occupied states. As a real space basis for the atomic orbitals </Wcr( r ) Slater-type orbitals [8] 
were used with effective charges taken from Ref. [9]. Table I lists the major nonvanishing Coulomb matrix elements 
which were used in this work. 


ELECTROSTATIC CONSIDERATIONS 

When depositing the junction on the metallic substrate coated by an insulating film, the parameters of the many- 
body Hamiltonian introduced in the previous sections are renormalized by image charge effects [10]. Also the presence of 
the other metallic lead (the tip) modifies the spectrum of the system. An additional source of complexity is introduced 
by the application of an external bias across the junction which brings the system out of equilibrium. Despite the 
complexity of a rigorous self-consistent treatment of the electrostatics of a nanojunction, the essential phenomena can 
be summarized to two: i) the image charges in the nearby metallic leads tend to stabilize an excess of charge on the 
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molecule, by reducing its addition energy; ii) The polarizability of the molecule implies that the bias voltage between 
source and drain is not entirely dropping at the contacts between the molecule and the leads. Depending on the 
geometry of the junction a considerable fraction can also drop across the molecule itself. Considering the weak coupling 
to the leads and the planar configuration of the CuPc in the STM set up, we have included the previous phenomena in 
the following minimal model. Firstly, we have introduced the effective Hamiltonian H mo i_ env = — 5{ C (N — No) 2 where 
N counts the number of (frontier) electrons on the molecule, No = 3 is the number of (frontier) electron associated 
to the neutral molecule and Ac is the strength of the image charge renormalization. Secondly, we have assumed the 
many-body levels obtained from the diagonalization of H mo i + Hmoi-env to be independent of the applied bias voltage. 
Moreover, for the bias dependance of the leads chemical potentials, we have used Rs/t = Ro ± e<ag / T Vb but requiring 
as + = 1- The values of Ac? and c^s/t/m are calculated according to the following electrostatic considerations. 

From the addition energy of the neutral molecule Uo = FW 0 +i,o — 27Av 0 ,o + En 0 -i,o we associate a capacitance to 
the molecule Cm = e 2 /Uo- For the tip-molecule and the substrate-molecule capacitances we adopt the parallel plate 
model and define Ct = e^A/h and Cs = eo e r A/d, where A = 144 A 2 is an estimate of the CuPc single molecule 
surface, h is the tip-molecule distance, e r = 5.9 is the relative permittivity of NaCl, d is the thickness of the NaCl 
thin film and eo is the vacuum permittivity. Connecting these three capacitances in series, we obtain an estimate of 
the relative potential drops 


^S/T/M — r (22) 

^S/T/M 

where C t “J = Cg -1 + C^ 1 + C^ 1 . The relative potantial drop on the molecule for h = 5 Aand d = 8.1 A(thickness 
of a trilayer NaCl) and Uo = 2.7 eV is about a quarter of the applied bias. 

The estimate for the image charge parameter Ac proceeds from the same model. First we calculate the electrostatic 
energy associated to the three capacitors Ct, Cm and Cs in series when: i) No external bias is applied but the first 
and the last plate are grounded and ii) a unit charge is deposited between Ct and Cm (Cm and Cs). We call the 
associated electrostatic energy E up (E^ own ). From a simple calculation one obtains: 


E, 


up 


2 + C T ’ 


E, 


down 


Cm+Cs 

Finally, for the image charge parameter we write: 


2 zife+Cs 


(23) 


- <5 ic = Eup + 2 Edown - U 0 /2, (24) 

i.e. the difference between the average electrostatic energy and the average energy needed to charge the isolated 
molecule. The average of the electrostatic energies gives an estimate of the energy needed to charge the molecule in 
presence of the leads. If we subtract from it the average addition energy of the isolated molecule Co/2, for which we 
already account in the many-body Hamiltonian (21), we obtain indeed an estimate of the image charge effects. 

Remarkably, using the model that we just exposed, we could fit the absolute spectroscopic position of the anionic 
and cationic transition and obtain, in accordance with the experiments, standard topographical images for CuPc on 
Cu[100] and trilayer NaCl [11], as well as CuPc on Cu[lll] and bilayer NaCl [12] within essentially the same set of 
fitting parameters e mo i and A for the isolated molecule, see Table I. 
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NaCl(3ML) NaCl(2ML) 


Cu[100] 

Cu[lll] 

00 (eV) 
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5.00 

d (A) 

8.1 

6.0 
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0.32 
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as 
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0.12 
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0.95 
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0.81 
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VT (V) 
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V£t (V) 

-2.72 
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